# Run regressions for survey analysis



# Setup ----
library(tidyverse)
library(lfe)
library(sandwich)
library(lmtest)
library(clubSandwich)
library(stargazer)

# Load survey data as "a"
# SOURCE: please directly contact the authors of You and Kousky (2024)






#####################################################################
###################### Regression Analysis ##########################
#####################################################################

## Propensity Score Weighting ------------------------------------------


match_glm <- glm(AnyIns ~ home_residents_num_new + 
                   log_home_tenure +
                   home_mortgage_dummy + 
                   renter +
                   single_family + 
                   non_children_disability_elderly_pets +
                   
                   hhi_low +
                   hhi_median +
                   employment_parttime_selfemployed +
                   employment_retired +
                   employment_Other +
                   nonwhite +
                   
                   damage_severity_home_num +      # home damage:            numeric 0-5 (0=no damage, 1=minimum damage..., 5=completely destroyed)
                   damage_severity_contents_num +  # home content damage:    numeric 0-5 (0=no damage, 1=minimum damage..., 5=completely destroyed)
                   evac_totaldollars_000 +         # evacuation costs:       $thousand
                   service_disrupt_cost_num +      # service disruption:     numeric 0-5 (0=no disrupt, 1=had disruption but no cost, 2-minor costs..., 5=Extreme Costs)
                   damage_car_Yes  +               # car damage:             0 or 1      (no car damage severity info except for michael)
                   lostincome_dummy +              # lost income:            0 or 1
                   costs_additional_Yes +          # additional costs outside evacuation:    0 or 1
                   prefer_not_answer +
                   
                   savings_ind +
                   Event.Name, data = a, family = "binomial")

a$predicted_AnyIns_glm = predict(match_glm, type="response")
a$weight_ATE_glm = ifelse(a$AnyIns==1, 1, a$predicted_AnyIns_glm/(1-a$predicted_AnyIns_glm))





#-------------------------------------------------------------------------------
# Define the control variables 
# Note that in the regressions we use "update()" to add the control variables 
#-------------------------------------------------------------------------------


# control for (1) disaster costs; (2) event dummy; (3) home characteristics; (4) household charactersitcis
controls  <- ~  damage_severity_home_num +      # home damage:            numeric 0-5 (0=no damage, 1=minimum damage..., 5=completely destroyed)
                damage_severity_contents_num +  # home content damage:    numeric 0-5 (0=no damage, 1=minimum damage..., 5=completely destroyed)
                evac_totaldollars_000 +         # evacuation costs:       $thousand
                service_disrupt_cost_num +      # service disruption:     numeric 0-5 (0=no disrupt, 1=had disruption but no cost, 2-minor costs..., 5=Extreme Costs)
                damage_car_Yes  +               # car damage:             0 or 1      (no car damage severity info except for michael)
                lostincome_dummy +              # lost income:            0 or 1
                costs_additional_Yes +          # additional costs outside evacuation:    0 or 1
  
                home_residents_num_new + 
                log_home_tenure +
                home_mortgage_dummy + 
                renter +
                single_family + 
                non_children_disability_elderly_pets +
  
                hhi_low +
                hhi_median +
                employment_parttime_selfemployed +
                employment_retired +
                employment_Other +
                nonwhite +
                
                prefer_not_answer









################################################################################
# Table 8: Effcts of Any Insurance
################################################################################


## unweighted ---------------------------------------------------------------

mylogit1 <- glm(update(controls, enough_money_low                ~ AnyIns + savings_ind + Event.Name + .), data = a,  family = "binomial")
mylogit2 <- glm(update(controls, financial_yr_post_veryworse_ind ~ AnyIns + savings_ind + Event.Name + .), data = a,  family = "binomial")
mylogit3 <- glm(update(controls, Unmet_Needs                     ~ AnyIns + savings_ind + Event.Name + .), data = a,  family = "binomial")




# Pseudo R^2 for each logitstic regressions 
Pseudo.Rsquared.1 <- with(summary(mylogit1), 1 - deviance/null.deviance)
Pseudo.Rsquared.2 <- with(summary(mylogit2), 1 - deviance/null.deviance)
Pseudo.Rsquared.3 <- with(summary(mylogit3), 1 - deviance/null.deviance)



# cluster by event
m1 = coeftest(mylogit1, vcov. = vcovCR(mylogit1, cluster = a$Event.Name, type = "CR2"))
m2 = coeftest(mylogit2, vcov. = vcovCR(mylogit2, cluster = a$Event.Name, type = "CR2"))
m3 = coeftest(mylogit3, vcov. = vcovCR(mylogit3, cluster = a$Event.Name, type = "CR2"))







## weighted --------------------------------------------------------------------
mylogit1_w2 <- glm(update(controls, enough_money_low                ~ AnyIns + savings_ind + Event.Name + .), data = a, weights = a$weight_ATE_glm,  family = "quasibinomial")
mylogit2_w2 <- glm(update(controls, financial_yr_post_veryworse_ind ~ AnyIns + savings_ind + Event.Name + .), data = a, weights = a$weight_ATE_glm,  family = "quasibinomial")
mylogit3_w2 <- glm(update(controls, Unmet_Needs                     ~ AnyIns + savings_ind + Event.Name + .), data = a, weights = a$weight_ATE_glm,  family = "quasibinomial")

Pseudo.Rsquared.1.w2 <- with(summary(mylogit1_w2), 1 - deviance/null.deviance)
Pseudo.Rsquared.2.w2 <- with(summary(mylogit2_w2), 1 - deviance/null.deviance)
Pseudo.Rsquared.3.w2 <- with(summary(mylogit3_w2), 1 - deviance/null.deviance)

m1_w = coeftest(mylogit1_w2, vcov. = vcovCR(mylogit1_w2, cluster = a$Event.Name, type = "CR2"))
m2_w = coeftest(mylogit2_w2, vcov. = vcovCR(mylogit2_w2, cluster = a$Event.Name, type = "CR2"))
m3_w = coeftest(mylogit3_w2, vcov. = vcovCR(mylogit3_w2, cluster = a$Event.Name, type = "CR2"))


Pseudo.Rsquared.survey <- c(Pseudo.Rsquared.1, Pseudo.Rsquared.1.w2, 
                            Pseudo.Rsquared.2, Pseudo.Rsquared.2.w2, 
                            Pseudo.Rsquared.3, Pseudo.Rsquared.3.w2)


stargazer( m1, m1_w, 
           m2, m2_w, 
           m3, m3_w, 
           
           
           type = 'text',
           #         type = 'latex',
           omit = c("Event.Name", "savings_ind", "num", "Yes", "dollars", "dummy",
                    "home_residents", "home_tenure", "home_mortgage", "renter", "single_family", 
                    "hhi", "employment", "nonwhite", "children", "prefer_not_answer"),
           omit.stat = c('ser', 'F'),
           omit.table.layout = 'n',
           add.lines = list(
             c( "Weighted", "No", "Yes", "No", "Yes", "No", "Yes"),
             c("Adj./Pseudo R-squared" , format(Pseudo.Rsquared.survey, digits = 3)
               )
           ),
           float = F,
           column.sep.width = "2pt",
           notes.append = F,
           font.size="small"
)







################################################################################
# Table 9: Alternative Funding Sources for Recovery
################################################################################

# Funding Sources ----------------------------------------
mylogit4_w2  <- glm(update(controls, fundingsource_Family_friends    ~ AnyIns + savings_ind + Event.Name + .), data = a,  family = "quasibinomial", weights = a$weight_ATE_glm)
mylogit5_w2  <- glm(update(controls, fundingsource_FEMA_grant        ~ AnyIns + savings_ind + Event.Name + .), data = a,  family = "quasibinomial", weights = a$weight_ATE_glm)
mylogit6_w2  <- glm(update(controls, fundingsource_SBA_loan          ~ AnyIns + savings_ind + Event.Name + .), data = a,  family = "quasibinomial", weights = a$weight_ATE_glm)
mylogit7_w2  <- glm(update(controls, fundingsource_Mysavings         ~ AnyIns + savings_ind + Event.Name + .), data = a,  family = "quasibinomial", weights = a$weight_ATE_glm)
mylogit8_w2  <- glm(update(controls, fundingsource_Myemployer        ~ AnyIns + savings_ind + Event.Name + .), data = a,  family = "quasibinomial", weights = a$weight_ATE_glm)
mylogit9_w2  <- glm(update(controls, fundingsource_charity_nonprofit ~ AnyIns + savings_ind + Event.Name + .), data = a,  family = "quasibinomial", weights = a$weight_ATE_glm)
mylogit10_w2 <- glm(update(controls, fundingsource_bank_loan         ~ AnyIns + savings_ind + Event.Name + .), data = a,  family = "quasibinomial", weights = a$weight_ATE_glm)
mylogit11_w2 <- glm(update(controls, fundingsource_credit_card       ~ AnyIns + savings_ind + Event.Name + .), data = a,  family = "quasibinomial", weights = a$weight_ATE_glm)
mylogit12_w2 <- glm(update(controls, fundingsource_Local_goverment   ~ AnyIns + savings_ind + Event.Name + .), data = a,  family = "quasibinomial", weights = a$weight_ATE_glm)
mylm13_w2    <-  lm(update(controls, fundingsource_number_new        ~ AnyIns + savings_ind + Event.Name + .), data = a,  weights = a$weight_ATE_glm)


n3  = coeftest(mylogit4_w2,  vcov. = vcovCR(mylogit4_w2,  cluster = a$Event.Name, type = "CR2"))
n4  = coeftest(mylogit5_w2,  vcov. = vcovCR(mylogit5_w2,  cluster = a$Event.Name, type = "CR2"))
n5  = coeftest(mylogit6_w2,  vcov. = vcovCR(mylogit6_w2,  cluster = a$Event.Name, type = "CR2"))
n6  = coeftest(mylogit7_w2,  vcov. = vcovCR(mylogit7_w2,  cluster = a$Event.Name, type = "CR2"))
n7  = coeftest(mylogit8_w2,  vcov. = vcovCR(mylogit8_w2,  cluster = a$Event.Name, type = "CR2"))
n8  = coeftest(mylogit9_w2,  vcov. = vcovCR(mylogit9_w2,  cluster = a$Event.Name, type = "CR2"))
n9  = coeftest(mylogit10_w2, vcov. = vcovCR(mylogit10_w2, cluster = a$Event.Name, type = "CR2"))
n10 = coeftest(mylogit11_w2, vcov. = vcovCR(mylogit11_w2, cluster = a$Event.Name, type = "CR2"))
n11 = coeftest(mylogit12_w2, vcov. = vcovCR(mylogit12_w2, cluster = a$Event.Name, type = "CR2"))
n12 = coeftest(mylm13_w2,    vcov. = vcovCR(mylm13_w2,    cluster = a$Event.Name, type = "CR2"))

Pseudo.Rsquared.4.w2  <- with(summary(mylogit4_w2),  1 - deviance/null.deviance)
Pseudo.Rsquared.5.w2  <- with(summary(mylogit5_w2),  1 - deviance/null.deviance)
Pseudo.Rsquared.6.w2  <- with(summary(mylogit6_w2),  1 - deviance/null.deviance)
Pseudo.Rsquared.7.w2  <- with(summary(mylogit7_w2),  1 - deviance/null.deviance)
Pseudo.Rsquared.8.w2  <- with(summary(mylogit8_w2),  1 - deviance/null.deviance)
Pseudo.Rsquared.9.w2  <- with(summary(mylogit9_w2),  1 - deviance/null.deviance)
Pseudo.Rsquared.10.w2 <- with(summary(mylogit10_w2), 1 - deviance/null.deviance)
Pseudo.Rsquared.11.w2 <- with(summary(mylogit11_w2), 1 - deviance/null.deviance)
Pseudo.Rsquared.12.w2 <- with(summary(mylogit12_w2), 1 - deviance/null.deviance)
Adj.Rsquared.lm.13.w2 <- summary(mylm13_w2)$adj.r.squared

Pseudo.Rsquared.source <- c(Pseudo.Rsquared.4.w2, Pseudo.Rsquared.5.w2,  Pseudo.Rsquared.6.w2, 
                            Pseudo.Rsquared.7.w2, Pseudo.Rsquared.8.w2,  Pseudo.Rsquared.9.w2, 
                            Pseudo.Rsquared.10.w2, Pseudo.Rsquared.11.w2, Pseudo.Rsquared.12.w2, 
                            Adj.Rsquared.lm.13.w2)

library(stargazer)
stargazer(n3, n4, n5, n6, n7, n8, n9, n10, n11, n12, 
          type = 'text',
          #         type = 'latex',
          column.labels = c("Family/Friends", "FEMA Grant", "SBA Loan", "Savings", "Employer", "Charity", "Bank Loan", "Credit Card", "Local Gov", "# of Sources"),
          omit = c("Event.Name", "savings_ind", "num", "totaldollars", "damage_car_Yes", "lostincome_dummy", 
                   "costs_additional_Yes", "home_residents", "home_tenure", "home_mortgage", 
                   "renter", "single_family", "hhi", "employment", "nonwhite", "children", "prefer_not_answer"),
          omit.stat = c('ser', 'F'),
          omit.table.layout = 'n',
          add.lines = list(
            c("Adj./Pseudo R-squared" , format(Pseudo.Rsquared.source, digits = 3))
          ),
          float = F,
          column.sep.width = "2pt",
          notes.append = F,
          font.size="small"
          
)









################################################################################
# Table 10: Perception of Insurance Usefulness and New Flood Insurance Purchases
################################################################################

r0  <- lm(update(controls,  help_ins_new ~ AnyIns + savings_ind + Event.Name + .),                                          data = a)
r1  <- lm(update(controls,  help_ins_new ~ savings_ind + Event.Name + .),                                                   data = a[a$AnyIns==1,])
r2  <- lm(update(controls,  help_ins_new ~ savings_ind + Event.Name + .),                                                   data = a[a$AnyIns==0,])
r3  <- glm(update(controls, ins_flood_switched ~ ins_home_rent_Yes + savings_ind + Event.Name + .), family = "binomial",    data = a[a$ins_flood_event=="No",])


r0_cluster = coeftest(r0, vcov. = vcovCR(r0, cluster = a$Event.Name, type = "CR2"))
r1_cluster = coeftest(r1, vcov. = vcovCR(r1, cluster = a[a$AnyIns==1,]$Event.Name, type = "CR2"))
r2_cluster = coeftest(r2, vcov. = vcovCR(r2, cluster = a[a$AnyIns==0,]$Event.Name, type = "CR2"))
r3_cluster = coeftest(r3, vcov. = vcovCR(r3, cluster = a[a$ins_flood_event=="No",]$Event.Name, type = "CR2"))


Adj.Rsquared.lm.r0 <- summary(r0)$adj.r.squared
Adj.Rsquared.lm.r1 <- summary(r1)$adj.r.squared
Adj.Rsquared.lm.r2 <- summary(r2)$adj.r.squared
Pseudo.Rsquared.r3 <- with(summary(r3), 1 - deviance/null.deviance)


Pseudo.Rsquared.perception <- c(Adj.Rsquared.lm.r0, Adj.Rsquared.lm.r1, Adj.Rsquared.lm.r2, 
                                Pseudo.Rsquared.r3)

library(stargazer)
stargazer(r0_cluster, r1_cluster, r2_cluster, r3_cluster, 
          type = 'text',
          #         type = 'latex',
          omit = c("missing", "Event", "savings_ind", "num", "log", "dummy", "costs", "dollar", 
                   "single_family", "renter", "house", "damage", "children", "employment", "nonwhite", "prefer_not_answer"),
          omit.stat = c('ser'),
          omit.table.layout = 'n',
          add.lines = list(
            c("Cluster by Event",    "Yes", "Yes", "Yes", "Yes"),
            c("Pseudo R-squared" , format(Pseudo.Rsquared.perception, digits = 2))
          ),
          float = F,
          column.sep.width = "2pt",
          notes.append = F,
          font.size="small"
          
)



























